* This do file replicates Nunn (2010) adding our controls to his specifications

* Set file directory path using <cd>


import excel "Figure 5\priests.xls", sheet("Survival data") firstrow case(lower) clear

destring death_yr_gha birth_year max_year min_year native, force replace

***** Bringing the data set in shape
// taking spells during which missionaries are not observed as "leaving the country"
* sort id max_year death_yr 

* drop if priestid==96 | priestid==133 | priestid==136 | priestid==2150 // Hayford, Hayfron and sons
drop if priestid==.

// Multiple visits - coded as different priests
replace priestid=priestid*1000000+_n if multiple==1 // Christian Protten, 3 times


// taking first/last year of observation as entry/exit into the Gold Coast
egen first_year=min(min_year), by(priestid)
egen last_year=max(max_year), by(priestid)
sort priestid first_year min_year
label var first_year "First year of service"
label var last_year "Last year of service"

** SVD Togo
replace death_yr_gha=death_yr_togo if church=="SVD Togo" & death_year<=1917 & death_yr_gha==.
replace last_year=1917 if last_year==. & church=="SVD Togo"

// Drop duplicates
duplicates drop priestid first_year last_year death_yr_gha, force 
drop min_year max_year
gen age_1st_station=first_year-birth_year
* reg age_1st_station first_year basel

** Coding died dummy variable indicating last year of survival in Ghana
* replace max_year=death_yr if death_yr!=.
replace last_year=death_yr_gha if last_year==.
gen died=(death_yr_gha==last_year)

* replace max_year=max_year+0.5 if min_year==max_year
replace last_year=last_year+0.5 if first_year==last_year

* Duration of stationing
gen duration=last_year-first_year
reg duration basel methodist catholic bremen if died==0 & native==0 & church!="SVD Togo"

* Recode multiple visits - coded as the same priests
replace priestid=3017 if priest=="Christian Protten"



********************************************************************************
* Declare data survival data
stset last_year, failure(died==1) origin(time first_year) id(priestid) time0(first_year)

* Generate time dummies to reflect the advance in medicine (quinine)
stsplit new, at(1(1)80)

* gen decade=0
* forval i=20(20)90 {
* replace decade=`i' if min_year>=18`i' & min_year<18`i'+20
*}
gen decade=.
forval i=20(10)90 {
replace decade=`i' if first_year>=18`i' & first_year<18`i'+20 & first_year>=1800
}

replace decade=80 if decade==90 // the 1890s decade only includes the year 1890, so we assign the year 1890 to the 1880s.

* Define Pre and post quinine era
gen postquinine=(first_year>=1850)
gen prequinine=(first_year<1850)
gen nonnative=(native==0)

replace died=0 if died==.
drop new

* gen age=min_year-birth_year
gen age=first_year-birth_year

gen age_missing=(birth_year==.)
replace age=25+_t0 if age_missing==1 // Average age at first station is 25

drop if first_year>1890
gen prequinine_native=prequinine*native
gen postquinine_native=postquinine*native
********************************************************************************


********************************************************************************
** Figure 5b *******************************************************************
********************************************************************************
gen age25=age-25

sts graph, by(postquinine native) adjustfor(age25) noshow tmax(10) ytitle(Proportion surviving) ytitle(, margin(medsmall)) title("") xtitle("Analysis Time, Years of service in Ghana") yla(.25(.25)1) xtitle(, margin(medsmall)) legend(order(1 "Pre-quinine x European" 2 "Pre-quinine x African" 3 "Post-quinine x European" 4 "Post-quinine x African") cols(2) ) 
* Line properties were changed with Stata's graph editor

graph export "Figure 5\Figure5b.pdf", replace

********************************************************************************
** Figure 5a ********************************************************************
********************************************************************************

** Number of European missionaries in Ghana (by denomination)
tabstat _st if native==0, statistics(count) by(first_year)

drop if native==1
collapse (sum) _st, by (first_year)
rename first_year year

rename _st n_europ
label var n_europ "N European Missionaries"

* Merge with death rate data
merge 1:1 year using "Figure 5\mortalityrates.dta", nogenerate // contains the death rate estimates

gen mr5=dr5/1000
label var mr5 "Mortality Rate"
twoway (lowess mr5 year if year<=1890 & year>=1750, lwidth(medium) xline(1850, lpattern(dash)) xlabel(1750(20)1890) ytitle("Mortality Rate")  yscale(range(0 0.4)) lcol(black)  lpattern(dash) lwidth(medthick)) (line n_europ year if year<=1890 & year>=1750, yaxis(2) ytitle("N European Missionaries", axis(2)) xtitle("") lcolor(black) legend(label(1 "Lowess Mortality Rate")  label(2 "N European Missionaries")))
********************************************************************************

* Add Pr-Quinine/Post-Quinine with Stata's graph editor

graph export "Figure 5\Figure5a.pdf", replace